require(clusrank)
require(readstata13)

d <- read.dta13("/Users/McKarlo/Dropbox/Research/Journal submissions/RoyalB/manuscript_data.dta")

# High- and low-outgroup exposure data frames
dh <- subset(d, village==3 | village==4 | village==5) 
dl <- subset(d, village==1 | village==2 | village==6) 
villageh <- dh[,1]
villagel <- dl[,1]

# Using list-subsetting to convert invidual columns into vectors

# Threat
t_commh <- dh[,6]
t_neighh <- dh[,7]
t_comml <- dl[,6]
t_neighl <- dl[,7]
clusWilcox.test(x = t_comml, y = t_neighl, cluster = villagel, 
                group = NULL, stratum = NULL,
                alternative = c("two.sided"),
                paired = TRUE, exact = TRUE, B = 2000, method = c("ds"))
clusWilcox.test(x = t_commh, y = t_neighh, cluster = villageh, 
                group = NULL, stratum = NULL,
                alternative = c("two.sided"),
                paired = TRUE, exact = TRUE, B = 2000, method = c("ds"))

# Cooperation
pgg_commh <- dh[,2]
pgg_neighh <- dh[,3]
pgg_comml <- dl[,2]
pgg_neighl <- dl[,3]
clusWilcox.test(x = pgg_comml, y = pgg_neighl, cluster = villagel, 
                group = NULL, stratum = NULL,
                alternative = c("two.sided"),
                paired = TRUE, exact = TRUE, B = 2000, method = c("ds"))
clusWilcox.test(x = pgg_commh, y = pgg_neighh, cluster = villageh, 
                group = NULL, stratum = NULL,
                alternative = c("two.sided"),
                paired = TRUE, exact = TRUE, B = 2000, method = c("ds"))


# Aggression
p_commh <- dh[,4]
p_neighh <- dh[,5]
p_comml <- dl[,4]
p_neighl <- dl[,5]
clusWilcox.test(x = p_comml, y = p_neighl, cluster = villagel, 
                group = NULL, stratum = NULL,
                alternative = c("two.sided"),
                paired = TRUE, exact = TRUE, B = 2000, method = c("ds"))
clusWilcox.test(x = p_commh, y = p_neighh, cluster = villageh, 
                group = NULL, stratum = NULL,
                alternative = c("two.sided"),
                paired = TRUE, exact = TRUE, B = 2000, method = c("ds"))


# Bias in aggressive individuals
# Aggressive and non-aggressive individuals data frames
da <- subset(d, aggressive==1)
dna <- subset(d, aggressive==0)
villagea <- da[,1]
villagena <- dna[,1]
pgg_comma <- da[,2]
pgg_neigha <- da[,3]
pgg_commna <- dna[,2]
pgg_neighna <- dna[,3]
clusWilcox.test(x = pgg_commna, y = pgg_neighna, cluster = villagena, 
                group = NULL, stratum = NULL,
                alternative = c("two.sided"),
                paired = TRUE, exact = TRUE, B = 2000, method = c("ds"))
clusWilcox.test(x = pgg_comma, y = pgg_neigha, cluster = villagea, 
                group = NULL, stratum = NULL,
                alternative = c("two.sided"),
                paired = TRUE, exact = TRUE, B = 2000, method = c("ds"))


